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Abstract 


A method is presented to calculate analytically the sensitivity derivatives of wing static 
aeroelastic characteristics with respect to wing shape parameters. The wing aerodynamic 
response under fixed total load is predicted with Weissinger’s L-method; its structural 
response is obtained with Giles’ equivalent plate method. The characteristics of interest 
in this study include the spanwise distribution of lift, trim angle of attack, rolling and 
pitching moments, induced drag, as well as the divergence dynamic pressure. The shape 
parameters considered are the wing area, aspect ratio, taper ratio, sweep angle, and 
tip twist angle. Results of sensitivity studies indicate that (1) approximations based 
on analytical sensitivity derivatives can be used for wide ranges of variations of the 
shape parameters considered, and (2) the analytical calculation of sensitivity derivatives is 
significantly less expensive than the conventional finite-difference alternative. 

Introduction 

During the design phase of an engineering system, numerous analyses are conducted to 
predict changes in the characteristics of the system due to changes in design variables. Usu- 
ally, this process entails perturbing each variable in turn, recalculating the characteristics, 
and evaluating the sensitivities with some sort of finite-difference process. The repeated 
analyses can drive the cost of design very high. An approach that has found increased 
interest recently in engineering design is analytical calculation of the sensitivity derivatives 
(Adelman and Haftka, 1987). Typically, the analytical approach requires less computa- 
tional resource than the finite-difference approach and is less subject to numerical errors 
(round-off or truncation). The analytical approach is best developed in parallel with the 
baseline analysis capability since it uses a significant portion of the numerical information 
generated during baseline analysis. 

In the design of modern aircraft, airframe flexibility is a concern from strength, control, 
and performance standpoints. To properly account for the aerodynamic and structural 
implications of flexibility, reliable aeroelastic sensitivity analysis is needed. Therefore, 
both structural and aerodynamic sensitivity analysis capabilities are necessary. 

Structural sensitivity analysis methodology has been available for well over two decades 
for both sizing (thickness, cross-section properties) and shape (configuration) variables 
(Adelman and Haftka, 1986). However, aerodynamic sensitivity analysis has been nonexis- 
tent until relatively recently. Some limited aerodynamic sensitivity analysis capability was 
developed for aircraft in subcritical compressible flow (Hawk and Bristow, 1984), but it only 
handled perturbations in the direction of the thickness of the wing (thickness, camber, or 
twist distribution). A new approach has been proposed by Yates (1987) that considers gen- 
eral geometry variations including planform for subsonic, sonic, and supersonic unsteady, 
nonplanar lifting-surface theory. 

Aeroelastic sensitivity analysis methodology has also been available for more than two 
decades for structural sizing variables. (See Haftka and Yates (1976) for an example 
application.) This is because changes in sizing variables affect exclusively the structural 
stiffness and mass distribution of the airframe and not its basic geometry. Therefore 
structural sensitivity analysis capability is sufficient. However, the lack of development in 
aerodynamic shape sensitivity analysis explains why there are very few results in aeroelastic 
shape sensitivity analysis. 

In a notable exception, Haftka et al. (1987) designed a sailplane wing under aeroelastic 
constraints and analyzed the design model with vortex lattice and finite element methods. 
A finite-difference aeroelastic sensitivity analysis capability is made possible by (1) devising 
a reduced order model to describe the wing static aeroelastic response and (2) using 
exact perturbation analysis to approximate changes in the vorticity vector with changes in 
the geometry. Because it retains techniques from both analytical and finite-difference 
sensitivity methodologies, this approach to sensitivity analysis should be described as 
semianalytical. 


The present study is a proof-of-concept that demonstrates the feasibility of calculating 
analytically the sensitivity of wing static aeroelastic characteristics to changes in wing 
shape. Of interest also is whether the curvature of the aeroelastic characteristics is small 
enough that analytical sensitivity derivatives can be used to approximate them without 
costly reanalyses. The flight regime is chosen to be subsonic and subcritical so that simple 
and inexpensive structural and aerodynamic analysis methodologies can be used. Yet the 
results are felt realistic enough that they may be used in conceptual design. It must be 
noted that the conclusions drawn here are strictly applicable to the analysis methodologies 
used here and may not hold true for other cases. 

This paper first describes the aeroelastic analysis methodology used. It combines 
Weissinger’s L-method (Weissinger, 1947) to predict the wing spanwise lift distribution 
with Giles’ equivalent plate analysis method (Giles, 1986) for structural analysis. The 
calculation of the sensitivity derivatives is described next. Finally, the methodology is 
used to investigate the sensitivity of a forward-swept wing to changes in wing area, aspect 
and taper ratios, and sweep and tip twist angles. The results are analyzed from the 
standpoints of both accuracy and computational cost. While the development of the 
sensitivity equations is given in the body of the text, the details of the derivation of the 
coefficients in those equations are given in appendices A-D. 

Nomenclature 


Ol 

chordwise position of the wing-box leading edge, normalized by the 
chord, measured from the leading edge 

a 2 

chordwise position of the wing-box trailing edge, normalized by the 
chord, measured from the leading edge 

A 

wing aspect ratio 

[A] 

aerodynamic matrix 

b 

wingspan 

Cr 

root chord 

ct 

tip chord 

{ccj} 

vector of product of section lift coefficients and chord length at aerody- 
namic control stations 

c la 

airfoil lift-curve slope 

[C] 

matrix as defined by equation (A10) 

CPU 

central processing unit 

d 

wing-box depth 

Di 

induced drag 

\D] 

matrix as defined by equation (A12) 

e 

position of airfoil center of pressure measured from the quarter-chord 
point, positive upstream, normalized by the chord 

M 

eigenvector, see equation (12) 

{*5,}, (4>> 

right and left eigenvectors at divergence 

E 

modulus of elasticity 

\E\ 

matrix as defined by equation (D2) 

{/} 

vector of applied point loads 
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G 

h(x,y ) 

M 

[K] 

L 

LE 

M 

M r 

n a 


n l 

n s 


N x 

Ny 

P 

q 

QD 

M 

s 

t 

TE 

T r 

M 

W] 

{w{x,y)} 

m 

[Wx] 


y 

Oio 

{a} 

{£} 

et 

V 


shear modulus 

transverse displacement of the wing surface 
identity matrix 
stiffness matrix 
lift 

leading edge 
flight Mach number 
rolling moment 

number of aerodynamic boundary-condition control points along 
semispan 

number of applied point loads along semispan 

number of terms in the approximation to the wing displacement 
function 

order of x in the polynomial displacement function h(x,y) 
order of y in the polynomial displacement function h(x, y) 
generic wing shape parameters 
dynamic pressure 
divergence dynamic pressure 

vector of amplitudes in the wing transverse displacement function 
(eq. (3)) 

wing area 

wing skin thickness 

trailing edge 

pitching moment 

vector of 1 

diagonal matrix of integration weights 
vector of displacement shape functions (eq. (3)) 

W^j is the value of entry j of vector {re} at the point of application of 
load i 

W x ij is the value of the streamwise derivative of entry j of vector {w} 
at control station i 

chordwise coordinate, positive downstream 
spanwise coordinate, positive out right wing 
wing root angle of attack 

vector of angles of attack at aerodynamic control stations 
vector of rigid twist angles at aerodynamic control stations 
wingtip twist angle, positive nose-up 

2y 

normalized spanwise coordinate, — 


3 



{0} vector of elastic twist angles at aerodynamic control stations 

A taper ratio 

A quarter-chord sweep angle, positive for sweepback 

v Poisson’s ratio 

Special Notations: 

{ } column vector 

[ ] square matrix 

( Y transposed quantity 

[ ] _1 inverse matrix 

( )' d{ )/dp 

Aeroelastic Analysis 

Aerodynamic Model 

The wing aerodynamic response is predicted by Weissinger’s L-method (Weissinger, 
1947), which was used herein as implemented for computations by DeYoung and Harper 
(1948). It is valid for moderate-to-high-aspect-ratio wings that are symmetric with respect 
to the root chord, have a straight quarter-chord line over each semispan, and have no 
discontinuities in twist. The airfoil section properties are assumed known and the flight 
regime may be compressible although it must be subcritical. In this method, the flow 
around the wing is modeled by a lifting line of vortices bound at the wing quarter-chord 
line. A no-penetration boundary condition is specified at n a control stations and that 
determines the spanwise distribution of vortex strength. Weissinger (1947) applies the 
boundary conditions at the three-quarter-chord point of each station; DeYoung and Harper 
(1948) modify that position to account for lift-curve slopes of less than the theoretical 2i r 
and also for effects of compressibility. A linear relationship between the vectors of local 
angles of attack and lift at the control stations results: 

{<*> = ( 1 ) 

The aerodynamic matrix [A] depends on the airfoil properties and the Mach number as 
well as the wing shape. Further details are given in appendix B. The total lift developed 
by the full-span wing is then given by: 

L = !*?{«}'[ V\{cc,} (2) 

Diagonal ( n a x n a ) matrix [V] contains shape- independent integration weights; its 
derivation is detailed in appendix A. 

Structural Model 

Giles’ equivalent plate analysis method assumes that the wing behaves like a plate and 
that its transverse displacements can be modeled by a polynomial in the chordwise and 
spanwise coordinates 

h{x,y) = {m(x,y)}*{s} (3) 

where vector {m} contains n s products of various powers of x and y. A careful selection 
of the exponents of the variables permits the specification of various types of boundary 


4 


conditions at the wing root. By applying the principle of virtual work, the wing static 
equilibrium equation .under a vector {/} of rq point loads is obtained 

MW = [wf {/} (4) 

The stiffness matrix [X] is consistent with the displacements given in equation (3); it is a 
function of the wing shape and sizing variables as well as wing material properties. This 
approach has shown very good results for both static and dynamic analysis of wings (Giles, 
1986 and 1987). It can handle fairly general planform geometries and boundary conditions, 
model complex wing cross-section geometries, include rib and spar caps, and permit the 
use of composite materials and the consideration of thermal loading. 

Aeroelastic Response 

To proceed with the calculation of the spanwise distribution of lift and the trim angle 
of attack under fixed lift, the vector of applied point loads is written in terms of the 
aerodynamic loads (see appendix A): 

{/} = jtiMW (5) 

The vector of angles of attack can be written as follows: 

{a} = a 0 {u} + {£} + {0} (6) 

All the entries of {u} are 1. Vector {0} can be derived from equation (3). At any control 
point i, the elastic twist is given by 


e i = -fa h i x vVi) 

For consistency with the aerodynamic model, the elastic twist is measured at the three- 
quarter chord. Then, 


W = ~[w x ]{s} (7) 

Combining equations (1), (4)— (7), as well as trim equation (2), we obtain the unknown 
angle of attack and the spanwise distribution of lift from: 

+ [V] —2b{u) 

§W‘M o 


Note that the inversion of the stiffness matrix precludes the use of free-body modes. 
Spanwise integration of the distribution of lift yields the rolling and pitching moments 



Mr = 

(9) 


Tr = [£>]{<*,} 

(10) 

The total induced drag is 

given by DeYoung and Harper (1948) as 



D, = ^W,}‘[£Hcc,} 

(11) 

Matrices [C] and [E] are independent on the wing shape; their derivation is given in 
appendices A and D, respectively. Matrix [D] is dependent on the wing shape; its derivation 
is given in appendix A. 
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One can obtain rigid-wing aerodynamic results for the analysis by setting q to zero in 
the left-hand side of equation (8) and performing the calculations of equations 9-11. 

Finally, the wing divergence dynamic pressure is found as the lowest eigenvalue q of the 
fixed-angle-of-attack problem: 


[A}+q 



[Wxwr'm'm 


{e} = {0} 


( 12 ) 


where {e} is the eigenvector. The matrices in equation (12) are nonsymmetric; therefore, 
the system has distinct right and left eigenvectors. 


Aeroelastic Sensitivity Analysis 

Sensitivity analysis begins with calculation of the derivatives of the distribution of lift 
and trim angle of attack. Taking the derivative of equation (8) with respect to p, we obtain 
after rearranging: 

" (-4 + f\WA\K\-'\wm 


— 2b{u} 
0 


ra' 


= f{2(6'{e} + 6{ £ }') j_ 


|/t]' + $ [f> 2 |W , x]|tf] _1 |W']‘] V] — 26'{u} 


l{u}*[V\ 


0 


/ w 1 1 


l 


a Q 


(13) 

As for all sensitivity equations, this equation is linear. The coefficient matrix on the left- 
hand side is identical to that of equation (8), which gives the distribution of lift and angle 
of attack. Further, it does not depend on the parameter with respect to which sensitivity is 
sought (a similar formulation is used by Yates (1987)). Therefore, only one calculation and 
inversion or factorization of the matrix is necessary for analysis and complete sensitivity 
analysis. There is one right-hand-side vector for each parameter of interest. The derivative 


of the matrix product [6^[W a; ][.K’] — *[W]*] is found by expansion. Calculations for the 
derivatives of the semispan b and the vector of twist {e}, the aerodynamic matrix [A], 
and the matrices [W x ] and [VF] are found in appendices A, B, and C, respectively. The 


derivatives of the flexibility matrix ([iF] *) are obtained from those of the stiffness matrix. 
Indeed, 


HC 1 = M 

After taking the derivatives and rearranging, 

[K)- v = -\K}- l {K]'[K\- 1 (14) 

In the present application, even though the sensitivity analysis is performed analytically, 
the derivatives of the stiffness matrix are found by finite difference. This practice makes the 
implementation almost as simple as for finite-difference sensitivity analysis while preserving 
much of the accuracy and cost advantages of the analytical approach. 

The shape derivatives of the rolling and pitching moments as well as the induced drag 
are obtained from equations (9)-(ll). We have 

K = |{u}'[V'|{2W/[C]{cc,} + 6 2 [C]{cc,}'} (15) 

Tl = | W'[V|{(/[Z>|{cc,} + HDJ'W + 6[C]{cc,}'} (16) 
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D'i = ^({cctfmcci) + {<*,}<[£]{«,}') (17) 

The details of the calculations for matrix [D]' are found in appendix A. Rigid-wing results 
can be obtained if the dynamic pressure is set to zero in equation (13). 

The derivatives of the divergence dynamic pressure are found from equation (12). If qj) 
is the divergence dynamic pressure and {e^,} is the right eigenvector, we have 


[ A \ + <?D 


r 6 2 




{ e £>} - M 


(18) 


Taking the derivative of this expression, premultiplying by the transposed left eigenvector 
, and rearranging, we find 


<4>i‘ 


[Vj 

{ e £>} 

{■ 

4>}‘ 


{ e D 

■} 


(19) 


Most of the matrix manipulations required in equation (19) were already performed while 
the derivatives of lift distribution and angle of attack were calculated. Therefore all that is 
necessary to evaluate the divergence dynamic pressure derivative is the calculation of the left- 
hand eigenvector. 


Applications 

The procedure described here was implemented on a Micro VAX II workstation in 
FORTRAN. Qualitatively, as can be inferred from the derivations, the development and im- 
plementation of the sensitivity analysis procedure are little more involved than those of the 
analysis procedure itself. This section of the paper describes some numerical applications. 

The wing model is described by five independent shape parameters: wing area S, aspect 
ratio A, taper ratio A, quarter-chord sweep angle A, and tip twist angle St ■ The model geometry 
is depicted in figure 1, and the problem parameters are given in table I. 

The number n Q of aerodynamic control stations was varied to examine convergence and 30 
was felt adequate; the effect of varying n a on accuracy and computational cost will be taken 
up in a subsequent section of this report. 

The wing structure is an aluminum wing-box of constant depth d and made of two flat cover 
skins of constant thickness f; it is located between the constant chord locations ai and «2- 

A measure of the flexibility of this wing is its divergence dynamic pressure of 16.3 kPa. 
For an air density of 1 kg/m 3 , the airspeed at which the elastic calculations are performed is 
89.4 m/sec, and the divergence speed is about twice as much at 180.3 m/sec. 

Results 


Sensitivity Results 

A sensitivity analysis is performed for the baseline design, and derivatives are generated 
analytically for rigid-wing and elastic-wing characteristics with respect to the five independent 
shape parameters. The characteristics include the rigid and elastic spanwise lift distribution 
(only wingtip values are given here, since they exhibit maximum aeroelastic effects), trim angle 
of attack, rolling and pitching moments, induced drag, and divergence dynamic pressure. These 
derivatives are used to approximate changes in the characteristics due to changes in the shape 
parameters. If r(p) is the value of characteristic r for a value p of a shape parameter, and if 
r'(p) is the sensitivity derivative of r at p with respect to p, then, for a small Ap, we can write: 
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r(p + A p) = r(p) + rf(p)Ap (20) 

In order to gauge the quality of the approximations, exact values of the perturbed characteristic 
r(p + A p) are generated by reanalysis for each value of A p. Figures 2-7 give examples of 
the sensitivity results so generated. All the figures share the same layout with the varied 
independent parameter on the horizontal axis and the characteristic on the vertical axis. 
Sensitivity results are given for both rigid- and elastic-wing cases. 

The results display two important properties. First, the straight lines of equation (20) are 
tangent to the corresponding sensitivity curves at the baseline values of the shape parameters, 
as, of course, they should be. Second, the curvatures of the characteristics are small enough 
so that sensitivity-derivative-based extrapolations can be used to approximate them for a wide 
range of variation of the shape parameters considered. The range of parameter variations for 
which the approximations are accurate enough depends largely on the application considered 
and the stage in the design process. Usually, however, the variations considered in a design 
effort are significantly smaller than those used in the present calculations, and the quality of 
the linear approximations should be quite adequate. Indeed, for parameter variations of ±10 
percent (±1° for the twist angle), the relative prediction error (((actual value — predicted 
value) /actual value |) never exceeds 1 percent in the results given here. Such a conclusion 
might differ in flight conditions that involve more complex flows as in the presence of shock 
waves or strong viscous interactions. 

Convergence of Calculated Parameters With Respect to n a 

The discretization of the aerodynamic model (n a ) is varied to assess its effect on convergence 
of the parameters and their derivatives. While n a is varied, the discretization of the structural 
model is held fixed (N x = 5, N y = 6) because, as indicated by Giles (1986), it cannot be 
increased significantly without risking singularities in the stiffness matrix. 

Table II shows the effect that aerodynamic discretization has on the convergence of the 
induced drag, the divergence dynamic pressure, as well as their derivatives. The induced drag 
and its derivatives converge more slowly than the divergence dynamic pressure. This is probably 
because the induced drag is very sensitive to variations in the shape of the load distribution 
and that the latter is slow to converge itself. For all the other characteristics considered, the 
trend is similar to that displayed by the divergence dynamic pressure. Note that, in both cases, 
the sensitivity derivatives of the characteristics converge about as fast as the characteristics 
themselves. 

Computational Cost 

Figures 8 and 9 show the effect that aerodynamic model discretization has on computational 
cost and its major components for generation of the sensitivity derivatives by finite difference 
and analytically. The cost is estimated in terms of both CPU time and total number of floating- 
point operations. The latter number is based on a double-precision LINPAK computing speed of 
.16 MFLOPS for a Micro VAX II workstation (Dongarra, 1985). For finite-difference sensitivity 
analysis, the cost only includes that of five full reanalyses (one for each independent shape 
parameter); no effort was made to avoid repeated calculations of invariant quantities. The 
cost of analytical calculation of the sensitivity derivatives is significantly lower than that of 
the finite-difference approach; the difference increases as the discretization is refined. For 
more elaborate analytical procedures that involve two-dimensional (surface panel) or three- 
dimensional (volume grid) discretizations, this effect is expected to become more pronounced. 

For both approaches, the cost associated with the calculations of the stiffness matrix and its 
derivatives remains constant since the structural model is unchanged. For the finite difference 
approach, the calculation of the aerodynamic matrix and its derivatives by finite difference 
is by far the major contributor to the total cost. For the analytical approach two items 
dominate the total cost for the finer discretization. First comes the cost of calculating the 
matrices [W] and [W x \, their derivatives, and the multiplications associated with changes 
between the discretization of the aerodynamic model and that of the structural model. Second 
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is that of the input /output operations relating to the transfer of information generated during 
analysis to the sensitivity analysis code. In fact the cost of computing the derivatives of the 
aerodynamic matrix is so small that it is not shown on the graph. 

Although the cost breakdowns discussed here may not be representative of those obtained 
with other analysis codes or other computer implementations, these results point to a cost 
advantage for the analytical approach to sensitivity analysis. 

Conclusions 

The results presented in this study show that, for subsonic subcritical flow, the curvatures 
of the wing aeroelastic characteristics are small enough that the characteristics can be well 
approximated by sensitivity-derivative-based linear extrapolations over ranges of variation 
of the shape parameters that are wide enough to be useful during the design process. 
Also, the analytical computation of the sensitivity derivatives is significantly less expensive 
computationally than the conventional finite-difference approach. In addition, with the simple 
aerodynamic and structural models used here, the derivation and implementation of the 
sensitivity analysis capability are little more complex than those of the corresponding analysis 
capability. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
March 29, 1988 
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Appendix A 


Aerodynamically Related Geometrical Parameters and Derivatives 

As shown in figure 10, the wing has a trapezoidal planform. The wing shape is 
completely described by five independent shape parameters: wing area S, aspect ratio A, 
taper ratio A, quarter-chord sweep angle A, and tip twist angle £f. The wing span b and 
the root and tip chords are dependent variables. They are found as 


b = \/~AS 

(Al) 

2 fs 

Cr= l + A VA 

(A2) 

2A [S 
C< = 1 + A V A 

(A3) 


The wing twist varies linearly between 0 at the wing root and at the wingtip. At a 
given spanwise location rj = 2 y/b, the twist angle is 

e(r/) = r)£ t (A4) 

Also, the wing chord c(rj) is given by 


c(h) = 


2^1 + (A - 1)??^ 
1 + A 



(A5) 


There are n a control stations in the aerodynamic model. They are defined so that their 
density increases toward the wingtip. The spanwise position of station i is 


2 y; 

Vi = ~r ~ = cos = cos 


VK 

2 n a 


(A6) 


Given the distribution {a} of angle of attack at those stations, the aerodynamic analysis 
determines the corresponding distribution of lift {cc z }. (See eq. (1).) For the sake of 
transferring the loads to the structural model, the distribution is replaced by a vector of 
point loads as illustrated in figure 11. At station i, the lifting load is given by 


fi = p cc li(Vi-i ~ *7t+l) 


(A7) 


with 


In vector form, we have 


fl = ^CC Z1 (1 -r/ 2 ) 
fn a = l^lnS'nria ~ X ) 


{/} = \bq[V){cci} 

with 


(A8) 
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[V] = diag [l - hi - h3> hi- 1 - hi+l> hn a “ l ] 



where [V] is a shape independent matrix. The total lift on both wings is obtained by 
summing up the point loads: 

Kcc,} (A9) 

where all the entries of {w} are 1. This equation is equivalent to the trapezoid integration 
rule. The moment of point load /,• with respect to the wing root is given by: 

b t 

m i = 7,mh 

so that the rolling moment is given by 


Afr = ^M*[V][C]{cc,} (A10) 

with 

[C] = diag [r)i,ri2, •••, Vi, •••, f?nj 

and matrix [C] is shape independent. Finally, as shown in figure 10, the pitching moment 
of fa with respect to the y-axis is (positive leading-edge up) 

fa x wifi (All) 

with 

x wi = J + Vi 2 tan A - ec(r? t -) 

where e , an airfoil property input into the calculations, is the distance between the center 
of pressure and the quarter-chord point normalized by the chord length; it is positive if the 
center of pressure is ahead of the quarter-chord point. The total pitching moment is then 
given by 


T r = ^{u}'(V][i>]{cc,> 

where 


(A12) 


\D\ — diag [ x w2i •••> x wi »•••> x wna\ 

The derivatives b ' , {e} 1 , and [£>]', with respect to the independent variables S , A , A, A, 
and St can be now traced through the expressions given in this appendix. 
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Appendix B 


Aerodynamic Matrix and Derivatives 

The derivation of [A] is detailed by DeYoung and Harper (1948). If the unknown 
distribution of circulation on the wing is approximated by n a terms of a Fourier series, 
specification of the no-penetration boundary conditions at n a control points (eq. (A6)) 
yields equation (1). [A] can be rewritten: 


[A) = \B)-±-\H)\L)\F] (Bl) 

The expressions in this appendix are given for symmetric loading; they include com- 
pressibility corrections and allow for input of section lift curve slope. Matrices [5] and [F] 
are shape independent; they are given in terms of the angles </> t - (eq. (A6)) by 


B ij — 2(6 'ij + bi2n a —j) 


= -2b, 


in n 


= 2b. 
= 2b 


n 


n a n a 

b = sin ^- ( i - (-iy 

(cos<f)j - cos<^) 2 l 4n a 


—i 


u _ n a 
**” 2sin& 


if * # j, j ^n a " 
i ± h j = n a 
i = j, j # n a 
i = j, j -n a > 

* 7 6 j 


(B2) 


and 


HI* 

ii 

if i 7^ T J ^ Wa 

— —f. . 

n a hj 

* = 1 , j + n a 

n a hj 

i # 1 , j = n a 

= 1 / • 
2 n a h > 

* = i, i = 

2n 0 -l 



fij = 2 H sin(*i0y) cos(f 1 </)j_ 1 ) 
i\odd 


(B3) 


Matrices [H] and [L] depend on the shape parameters exclusively through the following 
terms: 


Pi = 



(B4) 


r = tan = 


tan A 

\fi - M 2 


(B5) 
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where cj q and M are respectively the airfoil lift curve slope and the flight Mach number. 
Matrix [AT] is 


Matrix [A] is given by 


and 


Hj 


L tj 


:> •••) Pn a ] 

(B6) 

\ NS1 N32j 

(B7) 

) D22 { 


N1 ij = \J ( 1 + PiiVi ~Vj)rf + p^{rii - Vj ) 2 

N2 ij = \J (l + ?*(»?* - *lj ) T ) 2 + + Vj ) 2 

D22i = 1 + 2p i r H T 


D Uj = PiiVi - Vj) 


D2l ij = PiiVi + Vj) 


NS1 = 2 t 


W32 i = yj{ 1 + p 1 v 1 t) 2 + ~pfnf 


Because only matrices \H)' and [X] depend on the wing shape, the derivative of [A] with 
respect to any shape parameter p reads 

M' = -J^ | [»l'!il + [«1[£ |,|1 F1 (B8) 

The derivatives [H]' , can be determined from equations (Al), (A5), (B4), and (B6). For 
[A]', we have, for any shape parameter p 

= + /t3q\ 

dp dpi dp dr dp 

Derivatives dpjdp and dr/dp can be determined from equations (Al), (A5), (B4), and 
(B5). Now, if a stands for either or r, we have 


dh {j _ i dN ljj (Nlij - 1) dDlij i / N2jj \ dD2\ ij 

da m i;j da £>1?. da D2\\- \D22 l 1 ) da 

1 ( 1 dN2ij N2ij dD22 i \ N32 t dNSl 

+ D21ij \D22i da £>22? da j + D22 {j da 

N31 dN32j _ iV31iV32j dX722 ? - 
D22j da X>22? da 

1 i 
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and 


dNl 


dpi 


dNljj 


— = ((i + Pi(m - nj)r)(m - dj)r + - ^) 2 ) 

l ZJ 


dr 1 = TVT" ( ! + Pi{T,i ~ r? J )T ) Pi{Vi ~ ^ 


dD 1 


dpi 

dN 2 


— = (Vi - dj) 


dD 1 


<9r 


^=0 


= JV^ 0 1 + Pi(j1i ~ r] ^ T ^ r}l ~ *^ T + p ^ Vi + r? J^ 2 ) 
= (i + Piidi - dj)r) Piirii - dj) 


dN 2jj 

dr 

dD2l ij 

dpi 

dD 22j 

dpi 

dmij 

dpi 

dN 32.- 


“ZJ 


= (di + 

dj ) 

dD21ij 

dr 



dD22j 

= 2 diT 


i 

dr 

r\ 

dNSlj 

= u 


dr 

1 

~ JV32,- 

((1 + PidiT)d' t 

i T + Pidj 


= 0 


= 2 Pi r H 


dN 32, 1 fit x 

~dT = N32i {i + 


When evaluating terms on the diagonal of dL/dp j, care must be taken because the first 
term of equation (BIO) becomes indefinite for i = j. Using L’Hospital’s rule, it can be 
shown that 


if d i — dp then Jy[~ 


1 dNljj (Nljj - 1) dD 1 


ij dpj 


V 

DV- 


i] 




dpi 


0 
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Appendix C 


[W] and [Wx] Matrices and Derivatives 

The term W{j was defined previously as the value of the j th approximation function at the 
point of application of the ith load in the assumed displacement function of the equivalent plate 
analysis method. The j th approximation function reads as follows: 

wj = (— {JL-V yj (Cl) 

V^max / V 2/max / 

For wings clamped at the root, n x j takes all the values between 0 and N x ; for a given n x j , n y j 
varies between 2 and N y . The N x and N y values are chosen by the user; Giles (1986) suggests 
an upper limit of 7 for both of them. Also defined by the user are constant normalization 
lengths x max and 2/max- As detailed in appendix A, the aerodynamic loading is replaced by an 
equivalent set of point loads at the n a aerodynamic control stations. Therefore, we have for 
point i (see fig. 10) 


*wi = j + Vi 2 tan A - ec (*?i) 

6 

Vwi ~~ 2 ^ 


(C2) 


The parameters 6, c r , c(r/), and % are defined in terms of the wing independent shape variables 
in equations (Al), (A2), (A5), and (A6), respectively. If p is any of the independent shape 
parameters (5, A, A, and A) on which the structural model depends, then 


dW ij _ dw j fawi . dw j dywi fC3 s 

dp dx wi dp dy wi dp 

where dwj/dx w j and dwj/dy w i are calculated from equation (Cl). As to dx w i/dp and 
dy w Jdp, they are found by combining equations (C2), (Al), (A2), (A5), and (A6). This 
yields 


dx w i _ 1 - 4e(l - r?j) + (1 + A)r? t -Atan A -_4e?7;A 
~dS~ ~ 4v/l5(l + A) 

dxyi _ \/5(-l + 4e(l - rg) + 4er?tA) rg\/S tan A 
dA ~ 4^A 3 (1 + A) 4\/A 

dxwi _ Vjb 
dA 2 cos 2 A 

dxyjj = y/S(-2 + 8e(l-Vi)-8vi) 
d\ 4\/A(l + A) 2 


dy w j _ ms/A 
dS 4 y/S 

dywi _ mVs_ 
dA 4 y/A „ 

dyun q 

dA 

dywi j-. 

dX J 


(C4) 


The term W x ij was defined as the streamwise derivative of the y th approximation function at 
control station i, in the assumed displacement function of the equivalent plate analysis method. 
The streamwise derivative of the j’t h approximation function reads as follows: 


w jx = 


dWj 

n xj 

/ X \ n ^ 1 . 

( y \ n « 

dx 

^max 

V^max / 

\ 2/max/ 


(C5) 
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As mentioned in the discussion pertaining to equation (7), the elastic twist angle is measured 
at the three-quarter chord position of the aerodynamic model control stations. Therefore, we 
have for point i (see fig. 10) 


c W xi = t + Vi J tan A + Jc(%) 


Vwxi — Vwi “ 2 ^ 


(C6) 


The parameters b, c r , c(rj), and rjj are defined in terms of the wing independent shape variables 
in equations (Al), (A2), (A5), and (A6), respectively. If p is any of the independent shape 
parameters ( S , A, A, and A) on which the structural model depends, then 


dW, 


XIJ_ _ 


d dwj _ d 2 Wj dx n 


dp dp dx h 


dx 2 ■ 
ux wxi 


+ 


JPwj 


dy. 


WXl 


dp d Xxvii dy X wi 


dp 


(C7) 


d 2 Wj/dxl ]xi and d 2 Wj/dx wxi dy wxi are calculated from equation (C5). As to dx wx i/dp and 
dy wx i/dp, they are found by combining equations (C6), (Al), (A2), (A5), and (A6). This 
yields 


dx 


WXl 


ds 

d^wx 

dA 

U ' L WX l 

~dA 

U ' L "U) T 


d A 


l 3 — 2 Tji + 2f)i A + (l + A )t}iA tan A 

dVwxi 

dS 

_ ruy/A 
Ay/S 

4\/AS(l + A) 

i y/S(- 3-2(l-A)%) ! rn/S tan A 

dVwxi 

_ m/s 

4v^4 3 (1 + A) Ay/A 

d R 

Ay/A 

i _ Vib 
2 cos 2 A 

dy wxi 

dA 

= 0 

i y/S(-6 + 8r)i) 

dVwxi 

= 0 

4v / A(l + A) 2 

dA 


(C8) 
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Appendix D 
Drag Matrix 

DeYoung and Harper (1948) detail the calculations for the induced drag. In matrix form, 
the drag is given by 


D i= 


(Dl) 


Matrix [ E ] is independent of the wing shape; 

it is given 

equation (A6) and the coefficients b^j of equation (B2): 

Eij — 2 (bij + b{2n a ~j) s * n 

if *#/, 

= —2 bij sin (j>i 

if * # h 

= ~ipij + b&na-j ) sin 

if i^j, 

- 2ba sin <t>i 

if i = j , 

— bn a n a fina 

if i = j, 


in terms of the angles <j>i of 


* # n a , j ^ n a ' 
i # n a , j = n a 
i = n a > 

i 7^ n a 
i = n a 


(D2) 
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Table I. Baseline Model Parameters 


Wing geometry: 

S = 20 m 2 
A 1 = 7.5 
A = 0.5 
A = -20° 

£t= 0° 

Wing structural model: 

01 = 0.2 

02 = 0.7 

t = 0.02 m 
d = 0.1 m 

Airfoil properties: 

<H q =6 
e = 0 


Structural properties: 
E = 6.89 10* kPa 
v = 0.3 

G = 2.65 10 7 kPa 


Loading condition: 
q = 4 kPa 
M = 0.5 
L = 40000 N 

Analytical model discretizations: 
n a = ni= n s 
N x =5 
N y = 6 
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Table II. Effect of Aerodynamic Model Discretization on Convergence of Elastic 
Induced Drag and Divergence Dynamic Pressure and Their Derivatives 


Characteristics and 
derivatives 

n a = 10 

n a = 30 

n a = 50 

n 0 = 70 

N 

859.17 

852.70 

852.19 

852.05 

dD^/dS, N/m 2 .... 

-42.739 

-42.300 

-42.266 

-42.257 

dD?/dA, N 

-112.72 

-111.64 

-111.56 

-111.53 

dD\jd\, N 

24.489 

27.686 

27.936 

28.005 

dD?/dA, N/deg .... 

-0.10177 

-0.084786 

-0.084135 

-0.083973 

dD e Jde, N/deg .... 

2.3223 

2.8464 

2.8850 

2.8955 

<?£>, kPa 

16.308 

16.254 

16.250 

16.249 

dq D /dS, kPa/m 2 . . . 

-1.2220 

-1.2179 

-1.2176 

-1.2175 

dqp/dA, kPa 

-3.8220 

-3.8096 

-3.8086 

-3.8083 

dqo/dX, kPa 

-8.1577 

-8.1288 

-8.1265 

-8.1259 

dqp/d A, kPa/m 2 . . . 

6.7952 

6.7707 

6.7690 

6.7685 
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Figure 1. Baseline wing geometry. 
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Lift at tip 
section, 
N/m 



Figure 2. Sensitivity of lift at wingtip station to wing area. 


Angle of 
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Figure 3. Sensitivity of trim angle of attack to sweep angle. 
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Figure 4. Sensitivity of rolling moment to taper ratio. 



Figure 5. Sensitivity of induced drag to tip twist angle. 
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Figure 6. Sensitivity of induced drag to aspect ratio. 
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Figure 7. Sensitivity of divergence dynamic pressure to sweep angle. 
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Figure 8. Cost of finite-difference sensitivity analysis. 



Figure 9. Cost of analytical sensitivity analysis. 
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